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Abstract 

We study a discrete stochastic linear metapopulation model to understand the effect of risk 
spreading by dispersion. We calculate analytically the stable distribution of populations in different 
habitats. The simultaneous distribution of populations in habitats has a complicated self-similar 
structure, but the population in each habitat follows a log-normal distribution. A class of discrete 
stochastic matrix models were mostly dealt numerically. Our analytical predictions are robust in 
the wide range of parameters. Qualitative predictions of the current results should hold in the case 
of multiple habitats. We thus conclude that environmental stochasticity always promotes dispersal. 
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Environmental heterogeneity and flucuations are both important in ecological and evo- 



lutionary biology 



% Assuming that pop ul at 10 n S face thne-wying envi—al condi- 
tions, environmental heterogeneity often gives rise to paradoxical behavior [3J. For example, 
Jensen and Yoshimura [4| present a discrete-time model of offspring allocation into two 
habitats, either of which is so poor that populations cannot survive if the two habitats 
are exclusive. The authors demonstrate that dispersal can lead to population persistence. 
Similar paradoxical behaviors have been attractin g at tention in various research fields such 
as ecological biology js — 8 1 , financial economics js ill and information engineering 12]. In 



this paper, we study a discrete stochastic linear metapolulation model Q, which is given as 
discrete-time stochastic matrix models. This class of models have been mostly dealt numeri- 
cally by simulations, because of the analytical difficulty in tracking the stochastic processes. 
Here, we present an analytical method and show that these models are principally tractable 
analytically. To characterize long-term behavior of the entire system, geometric mean of 
local growth rate has been often used 0, Q, [ll|. However, this simple approach is not valid 
except for well-mixed cases. We present an accurate method to assess the long-term growth 
in this model. We derive analytically the stable distribution of populations. 

Let us consider populations that inhabits in n discrete habitats. Let Xi(t) be the number 
of individuals in patch % at time t. Thus, the state of the populations is described by a 
vector Xj(t) = (xi(t),x 2 (t), . . . x„(t)) T , where the superscript T represents the transpose. In 
each habitat, the population reproduces at random growth rates. Then a fraction of the 
population disperses from a habitat to another habitat. The population dynamics is given 
by a discrete-time stochastic matrix model 

x(i+l) = M(t)x(i), (1) 

where D is a time-independent dispersal matrix, and M(t) is a time-dependent diagonal 
matrix of local growth rate. In this paper, we focus on the case of two habitats for simplicity, 
but the results can be extended to a general case (i.e., n > 2). The dispersal matrix is given 

as 

n={ l ~ q qs ). (2) 

y qs 1 - q J 

Here the parameters q and s are between and 1. The migration rate q represents the 
proportion of population that migrates from a habitat to the other habitat. When q = 



0, the two habitats are isolated completely. The parameter s is survival rate during the 
transportation between habitats. The matrix of growth rate is written as 

MWH miW " I- (3) 
m 2 (t) 

We assume that the local growth rate mj(i) (i — 1 or 2) is a stochastic variable which 
takes one of two values, m_ with probability p and m + with probability 1 — p. Here, we set 
m_ < 771+. This stochastic fluctuation of the local growth rate comes from environmental 
fluctuations. Here, we neglect the temporal correlation of the local growth rates. But we 
take into account the correlation between two habitats. We denote by c the correlation 
coefficient between m\{t) and 7772 (i) at same time. If we denote by p ++ the probability that 
mi(t) = m + and 7772 (£) = rn + , we obtain 

P++ = (1 -p) 2 + cp(l -p) 
P— = P 2 + cp(l -p) 

(4) 

p+- = (1 - cM 1 -p) 

P-+ = (1 - c )p(! -P)- 

The initial populations at t = is set as (x\(0), X2(0)) T = (1, 1) T . 

To solve eq. (Q, we consider the dynamics of the ratio of populations in two habitats 



The dynamics of r(i) is written as 

mx(rj) 



Kt + l) = /f^ir(t)V (6) 



where /(r) is defined as 



/(r) = r(1 " (7) 



Equation (jHJ) indicates that the future value r(t + 1) depends only on the present value 
r(i), that is, r(t) follows one-dimensional Markov process. Thus, statistical state of r(t) is 
characterized by its density distribution ptij). The evolution of ptir) is described by the 
Frobenius- Perron equation [kJ Q. Using the probability p , p ++ defined in eq. (j4j), 



we have 

Pt+i(r) = P+-—g{r)pt(—f~ 1 {r)) 



m + \ m + ) 

+ (p++ +P— )g(r)pt(f~ 1 (r)), 

where g(r) is the derivative function of / _1 (r): 

l-2g-g 2 (s 2 -l) 

#( r ) = -7 7 7Y2 — • ( 9 ) 

(rqs + q — l) z 

The density distribution pt(r) converges to the stationary distribution p*(r) (or invariant 
measure) eventually, regardless of where it begins. The distribution p*(r) is determined 
by the recursive formula obtained by substituting p*(r) for pt(r) and p t+ i(r) in eq. (JHJ). 
Figure 1(a) shows an example of p*(r) with a logarithmic scale. This distribution p*(r) 
has a complicated self-similar structure. Figure 1(b) shows the time distribution obtained 
by simulation over 100000 time steps for the same parameter of Fig. 1(a). These two 
distributions agree exactly. Because of the ergodicity of a Markov process, the time average 
is equal to the ensemble average. The distribution of lnr(t) is symmetric around the vertical 
axis lnr(t) = (Fig. 1), because this model is invariant under the permutation of r(t) and 
l/r(t). This means that the distribution of r(t) coincides with that of l/r(t). In the special 
case qs = 1 — q, r(t) becomes one deterministically for t > 0, that is, p*(r) = S(r — 1). This 
case represents well-mixed populations, which were examined by prior theoretical works 

Q,flQ,Q. 

The dynamics of the populations x\{t) and ^(t) is described as 

x x {t + 1) = s 1 (0K(t)(l -q) + m 2 {t)qsr(t)}, (10) 
x 2 {t + 1) = x 2 (t)[m 2 (t)(l -q) + mi(t)qs/r(t)}. (11) 

A very important point is that the stochastic variable r(t) is determined by eq. (jHJ) and can 
be treated independent of xi(t) and x 2 (t). Thus, eqs. ( |T0|) and ( ITTjl are regarded as random 
multiplicative processes for xi(t) and x 2 (t), respectively. Taking into account the symmetry 
of r(t) and l/r(t), it is obvious that the two processes (fTUI) and (iTTj) are identical. Hence, 
it is sufficient to consider only ffTOl) . Taking the logarithms of both sides of eq. ( TTUT) and 
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FIG. 1. Stationary distribution of r(t) = X2{t)/x\{t) for m+ = 3, m_ = 0.01, c = 0, p = 0.3, 
q = 0.2, and s = 0.5. The horizontal axis is a logarithmic scale, (a) Theoretical result obtained by 
solving the Frobenius- Perron equation ([8]) numerically, (b) Time distribution obtained by computer 
simulation over 100000 time steps after a transient of 1000 time steps. We observe a good agreement 
between these two distributions. 

summing them from t = to t = T — 1, we obtain 

t=T-l 

lnxx(T) = lnxi(t + 1) — lnxi(i) 

dT°-i (12) 
= ln[mi(t)(l - q) + m 2 (t)qsr(t)\. 

t=o 

Thus, lnxi(T) is given by the sum of the time series of the effective growth rate ln[mi(i)(l — 
q) + m 2 (t)gsr(t)]. Because the effective growth rate depends on r(t), the effective growth 
rates have a temporal correlation. However, the auto-correlation function of the effective 
growth rates decays rapidly, as shown in Fig. 2(c). Thus, applying the central limit theo- 
rem, we expect the distribution of lnxi(T) after a long time (T ^> 1) approaches a normal 
distribution. Thus, X\(T) follows a log- normal distribution. The mean value of lnxi(T) 
is calculated by using the stationary distribution of r(t). Moreover, in an approximation 
in which the temporal correlation of r(t) is neglected, we can also estimate the variance of 
\nxi(T). In Figs. 2(a) and (b), we show examples of the evolution of the actual ensemble 
distribution and the approximation results. The averages (or peaks) of the actual distribu- 
tions agree exactly with the theoretical results. The actual variances are also equal to the 
theoretical results (Fig. 2(b)) when the auto-correlations are negligible (Fig. 2(c)). How- 
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FIG. 2. Evolution of probability distribution of \nx\(t) for (a) q = 0.2 and (b) q = 0.65, when 
t = 1000,2000, and 4000. The other parameters are m+ = 3, m__ = 0.01, c = 0, p = 0.3, and 
s = 0.5. The crosses stand for the distribution obtained by 100000 stochastic realizations. The 
lines stand for the approximation with neglecting time correlation of the effective growth rates, 
(c) The temporal correlation functions of the effective growth rates of x\{t) in the cases (a) and 
(b). We find a slight negative correlation at short time lags and rapid decay to zero. The negative 
correlation implies that the variance of the actual probability distribution tends to be smaller than 
the approximation with neglecting temporal correlation of the effective growth rate. (Color online) 



ever, the actual variances are smaller than the approximation results (Fig. 2(a)), when the 
auto-correlations are not neglected (Fig. 2(c)). Here the deviations are universally smaller 
because the growth rates tend to have a negative correlation in short-range intervals of time. 

Because X\{t) and X2{t) obey the same stochastic process, lna^CO follows the same 
normal distribution as \nxi(T). At a glance, the difference lna^CO — lnxi(T) looks like 
to follow a normal distribution. But this insight is not correct. Recall that In r(T) = 
ln^2(T) — lnxi(T). We have concluded that this distribution has a complicated form (as is 
seen in Fig. 1). This indicates that the two variables X\(T) and X2{T) are not independent of 
each other but have a complex relationship. Thus, the simultaneous distribution of X\(T) and 
Xi(T) has a complicated self-similar structure. For the well-mixed case (qs = l—q), xi(t) and 
X2(t) simply coincide. However, unless the populations are well-mixed, lnxi(T) +lna;2(T) 
does not follow a normal distribution and furthermore resulting Xi(T) + £2 CO does not 
follow a log-normal distribution. The current result holds regardless of the shape of the 
distribution of mi(t) and m2(t), since it is not because of the restriction that mi(t) and 
m,2(t) can have only two values. 
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Which population survives is assessed by the ensemble average of the long-term popula- 
tion growth 



For the well- mixed case (qs = 1 — q), the long-term growth rate (1131) is rewritten by the 
average of local growth rates 



Contrarily, for the isolated case (q = 0), it is calculated as (1 — p)lnm + +plnm_. For 
general cases, however, the long-term growth rate ffl3|) cannot be written in a simple form. 
To obtain the value of (113j) . we need to calculate (jl2jl by the stationary distribution of (JE]) 
with the help of computer. 

Figure 3 gives a comparison between numerical simulations and analytical results. If 
the long-term growth rate (fl3|) for q = is negative, the population cannot survive in a 
single habitat. As shown in Fig. 3(a), the growth rate increases with the migration rate 
q until the optimal migration rate q*, where the long-term growth rate has the maximal 
value. Consequently, the dispersal is advantageous for the populations to persist. Figure 
3(a) shows that the effect of dispersal is strong (weak) when the environments of habitats 
have a negative (positive) correlation. In Fig. 3(b), we plot the optimal migration rate q* 
as a function of survival rate s for five values of p. The curves for p and 1 — p coincide, 
because of the symmetry between the two habitats. Although the optimal migration rate 
q* decreases with decreasing s, q* remains finite even for relatively small s. Note that the 
average (4 lnxi(T)) is the logarithm of the geometric average of Xi(T). Since Xi(T) follows 
a lognormal distribution, the geometric average of X\(T) coincides with the median of X\(T). 
This means that when the long term population growth (eq. ffl3|) ) is positive (negative), the 
population grows (decays) with a probability of more than half. 

Risk-spreading phenomena have been treated in biology and economics applying dis- 
crete stochastic linear metapolulation models. Most these studies are based on numerical 
approaches, except well-mixed cases that can be solved by simple algebra. Here we have 
solved analytically the cases of non-well-mixed populations. Our approaches can be widely 
applicable for forecasting long term trends not only in physical phenomena, but also in 
population biology and economic forecasting. Even complex stochastic problems may be 




(13) 





(14) 



+2p + _ ln(m + + mJ) + ln(l — q). 
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(a) 



(b) 




FIG. 3. (a) Long term population growth rate is shown as a function of migration rate q for three 
values of correlation (c = —0.2,0,0.2). The other parameters are set as m + = 3, m_ = 0.01, 
p = 0.3, and s = 0.5. (b) The optimal migration rate is shown as a function of survival rate 
through migration s for p = 0.2, 03, 0.5, 0.7, 0.8. The curve is obtained by (|12p and the stationary 
distribution p*{r). The crosses represent the numerical result of (|13p for T = 10000, averaged over 
10000 stochastic realizations. 

tractable and solved analytically in this approach. 

We analyzed the case of two habitats for simplicity. In the case of more than two habitats, 
the qualitative prediction should be same as in the case two habitats. Here we could expect 
the same The above results are expanded to the case of networks with more than two 
habitats. Here, the ratios of populations among habitats have a complicated self-similar 
stationary distribution, while each population always has a log-normal distribution. In 
other words, the simultaneous distribution has a complicated structure, but its marginal 
distributions always follow a log-normal distribution. However, the numerical calculations 
become terribly cumbersome even in the case of three habitats. We should note here that this 



problem of multiple habitats may be treated in the framewor 



because some analytical solutions are already acquired 



13, 



£ of metapopulation dynamics 



181 ] . We may also apply the 



current model to the studies of horizontal gene transfer in two distinctive environmental 



states 



19 



201 ] . We can also add a nonlinear effect to the current linear model. Then, 
the stable distribution departs from a log-normal distribution. For example, the density 
dependent effect prevents a population from increasing to an infinitely large number. On 
the other hand, it is well-known that a finite injection leads to a power law distribution 
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21 



23j . Although these cases have not been investigated here, we can conjecture that the 



simultaneous distribution may still show have complicated self-similarity. In addition, some 
environmental factors are naturally correlated in successive times, resulting in the temporal 
correlations in local growth rates. These cases are treated numerically |7j], because temporal 
correlations are extremely complex and highly tedious to treat analytically. The current 
method may be applicable to solve these models analytically and we expect no qualitative 
differences from the current results. We also find robustness in environmental parameters 
m + , m_, p and c. We generally could conclude that the stochasticity always promotes 
dispersal to ensure the long term survival even if the cost of migration is considerably high. 
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